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Abstract 

We  propose  a  novel  method  to  implicitly  two-way  couple  Eulerian  compressible  flow  to  volumetric  Lagrangian 
solids.  The  method  works  for  both  deformable  and  rigid  solids  and  for  arbitrary  equations  of  state.  The 
method  exploits  the  formulation  of  [11]  which  solves  compressible  fluid  in  a  semi-implicit  manner,  solving 
for  the  advection  part  explicitly  and  then  correcting  the  intermediate  state  to  time  using  an  implicit 
pressure,  obtained  by  solving  a  modified  Poisson  system.  Similar  to  previous  fluid-structure  interaction 
methods,  we  apply  pressure  forces  to  the  solid  and  enforce  a  velocity  boundary  condition  on  the  fluid  in 
order  to  satisfy  a  no-slip  constraint.  Unlike  previous  methods,  however,  we  apply  these  coupled  interactions 
implicitly  by  adding  the  constraint  to  the  pressure  system  and  combining  it  with  any  implicit  solid  forces 
in  order  to  obtain  a  strongly  coupled,  symmetric  indefinite  system  (similar  to  [17],  which  only  handles 
incompressible  flow).  We  also  show  that,  under  a  few  reasonable  assumptions,  this  system  can  be  made 
symmetric  positive-definite  by  following  the  methodology  of  [16].  Because  our  method  handles  the  fluid- 
structure  interactions  implicitly,  we  avoid  introducing  any  new  time  step  restrictions  and  obtain  stable  results 
even  for  high  density-to-mass  ratios,  where  explicit  methods  struggle  or  fail.  We  exactly  conserve  momentum 
and  kinetic  energy  (thermal  fluid-structure  interactions  are  not  considered)  at  the  fluid-structure  interface, 
and  hence  naturally  handle  highly  non-linear  phenomenon  such  as  shocks,  contacts  and  rarefactions. 


1.  Introduction 

Direct  numerical  simulations  (DNS)  are  often  used  to  study  the  interactions  between  fluid  flows  and 
solid  structural  models.  Under  certain  assumptions  these  can  be  reduced  to  a  one-way  coupled  system;  for 
example  if  one  wishes  to  determine  the  steady-state  lift  of  an  airfoil  in  subsonic  flow,  it  is  often  reasonable 
to  simulate  the  airfoil  as  a  kinematic  body.  With  a  clever  choice  of  boundary  conditions,  one  can  even 
begin  to  examine  two-way  coupled  interactions,  albeit  in  a  limited  fashion.  In  the  more  general  case,  these 
assumptions  miss  the  interesting  two-way  coupled  interactions  between  the  fluid  and  the  structure.  These 
two-way  coupled  interactions  can  be  quite  important  and,  if  not  properly  captured  in  the  DNS,  can  lead  to 
non-physical  results.  It  is  therefore  important  to  have  a  robust  numerical  method  that  accurately  captures 
two-way  coupled  interactions  across  a  fluid-structure  interface. 

Methods  to  capture  fluid- structure  interactions  can  be  broadly  separated  into  two  categories.  Weakly 
coupled  (partitioned)  systems  interleave  the  disparate  subsystems  by  integrating  them  forward  in  time  sepa¬ 
rately,  using  each  others’  results  as  boundary  conditions  in  an  alternating  one-way  coupled  fashion  (see  e.g. 
[21,  15,  6]).  This  approach  is  appealing  as  it  permits  the  use  of  specialized  numerical  methods  for  each  of  the 
different  materials  with  only  slight  modifications  to  account  for  the  modified  time  integration  and  changing 
boundaries.  There  are  disadvantages  to  this  approach,  however,  for  example  new  and  poorly  understood 
stability  restrictions  arise  independent  of  the  individual  subsystems,  such  as  the  lumped- mass  instability 
discussed  in  [4].  The  alternative  is  to  employ  a  strongly  coupled  (monolithic)  system,  which  are  systems 
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where  the  fluid  and  structure  are  evolved  forward  in  time  simultaneously  using  a  solver  specially  crafted  to 
incorporate  phenomena  from  both  fluid  and  solid  phases.  Our  method  is  a  hybrid  of  the  two;  the  explicit 
components  of  both  fluid  and  solid  solvers  are  evolved  forward  independently,  while  the  implicit  components 
and  interactions  are  coupled  together  in  a  monolithic  solve. 

State-of-the-art  solvers  typically  use  an  Eulerian  framework  to  treat  fluid  flows  and  a  Lagrangian  frame¬ 
work  to  treat  solids,  and  so  any  coupled  system  must  do  one  of  three  things:  model  the  solid  in  an  Eulerian 
framework,  model  the  fluid  in  a  Lagrangian  framework,  or  find  a  way  to  couple  Eulerian  fluids  with  La¬ 
grangian  solids.  The  first  two  options  are  undesirable  as  they  impose  significant  limitations  on  the  numerical 
method,  for  example  Eulerian  models  only  capture  material  properties  (rather  than  tracking  them)  which 
makes  it  difficult  to  compute  time  history  variables  important  to  structural  simulation,  such  as  loading  and 
damage.  Many  fluid  Lagrangian  models  have  difficulty  in  obtaining  the  correct  shock  speeds  due  to  the 
lack  of  discrete  flux  differencing,  and  therefore  resort  to  artificial  viscosity  methods  that  require  a  number 
of  zones  within  a  shock  in  order  to  obtain  the  right  speed  [2,  3].  Lagrangian  fluid  models  also  struggle 
with  high-speed  and  deforming  flows,  as  large  deformations  can  cause  significant  numerical  errors  in  the 
flow  field  and  can  drive  the  time  step  to  zero.  This  can  be  partially  alleviated  by  applying  complex  and 
expensive  remeshing,  but  if  the  flow  field  tangles  and  inverts,  the  simulation  can  cease  altogether.  Arbitrary 
Lagrange-Eulerian  (ALE)  methods  address  the  problem  of  a  deforming  Lagrangian  fluid  grid  by  permitting 
the  fluid  grid  to  move  at  some  velocity  other  than  the  velocity  of  the  fluid,  but  this  can  still  lead  to  high 
aspect  ratios  that  necessitate  remeshing,  especially  in  the  presence  of  a  fluid-structure  interface.  We  address 
the  challenge  of  coupling  Eulerian  fluids  with  Lagrangian  solids  by  introducing  an  interpolation  operator, 
which  conservatively  maps  quantities  from  Eulerian  boundaries  to  nearby  Lagrangian  boundary  nodes,  and 
vice  versa. 

At  the  fluid-structure  interface  there  is  a  transfer  of  information.  This  information  transfer  can  be  handled 
by  weakly  coupling  each  separate  subsystem  using  a  one-sided  estimate  of  the  transfer,  or  by  strongly  coupling 
subsystems  together  and  introducing  new  variables  to  the  equations.  Weakly  coupled  approaches  have  been 
shown  to  give  high-fidelity  results  [1,  8,  7],  but  can  struggle  when  applied  to  a  system  with  high  density-to- 
mass  ratios  (and  are  prone  to  going  unstable,  as  we  discuss  in  Section  4.3).  These  problems  can  be  alleviated 
by  using  a  better  estimate  of  values  at  the  interface,  as  suggested  by  [12],  but  this  typically  involves  solving 
expensive  general  Riemann  problems  at  every  fluid-structure  face.  These  problems  can  be  avoided  entirely 
by  handling  the  interface  in  a  strongly  coupled  fashion,  but  previous  work  has  been  limited  to  incompressible 
flows  [14,  17].  Our  method  exploits  the  structure  of  [11],  which  treats  the  pressure  flux  of  compressible  flows 
implicitly.  This  permits  us  to  treat  the  fluid  pressure  as  an  implicit  force  on  the  solid,  and  use  an  implicit 
velocity  boundary  condition  on  the  Poisson  solve,  much  like  previous  strongly-coupled  work. 

Our  fluid  evolution  is  comprised  of  two  steps:  an  advection  stage  and  a  pressure  solver  phase.  This  permits 
us  to  address  the  complexities  arising  from  the  truly  non-linear  components  of  the  flow  separately  from  the 
linearly  degenerate  components.  In  the  pressure  phase,  we  freeze  everything  to  their  time  location  and 
perform  an  implicit  solve  for  the  fluid  pressure  and  solid  velocity.  It  is  in  this  phase  that  we  handle  the 
transfer  of  momentum  and  kinetic  energy  across  the  fluid-structure  interface,  and  as  such  it  is  important 
to  be  conservative  in  transferring  information  between  the  two  sets  of  degrees  of  freedom.  In  the  advection 
stage  no  information  should  be  transmitted  across  the  interface,  but  instead  we  must  address  the  issues  which 
arise  by  virtue  of  a  moving  solid  (i.e.  the  covering  and  uncovering  of  fluid  cells).  There  are  many  examples 
of  how  to  address  these  problems  in  the  literature,  for  example  we  could  track  cut  cells,  re-discretize  the 
fluid  in  an  ALE  formulation — all  of  which  significantly  complicate  the  fluid  evolution.  Instead  we  make  the 
key  observation  that  since  the  interface  is  a  contact  discontinuity  we  can  afford  to  be  non- conservative,  but 
only  in  the  linearly  degenerate  components  of  the  flow. 

In  a  traditional  explicit  method  the  linearly  degenerate  and  truly  non-linear  fluxes  aren’t  separated, 
and  as  such  these  methods  need  to  deal  with  all  of  the  complexities  of  moving  boundaries  and  information 
transmission  at  the  same  time.  That  is,  they  need  to  be  conservative  when  dealing  with  information  that 
crosses  the  interface  while  at  the  same  time  dealing  with  an  interface  that  moves.  Einally,  the  flux  needs  to  be 
re-examined  carefully  in  order  to  determine  what  forces  should  be  applied  to  the  interface.  One  could  modify 
traditional  methods  by  separating  the  conserved  quantities  into  their  Riemann  invariants,  and  be  conservative 
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in  the  truly  non-linear  invariants  while  allowing  the  linearly  degenerate  invariants  to  be  non-conservative 
— however  this  doesn’t  address  the  moving  boundary,  and  still  leaves  us  with  the  (poorly-understood)  CFL 
restriction  that  arises  from  explicit  fluid-structure  interactions.  Because  of  these  complications,  our  method 
hinges  on  the  existence  of  [11]. 


2.  Semi- implicit  compressible  flow 

We  briefly  describe  the  semi-implicit  evolution  for  compressible  flow  [11]  which  forms  the  basis  for  our 
implicit  coupling  scheme.  Consider  the  multi-dimensional  Navier-Stokes  equations,  given  by: 

(p\  (  V.pn  \  /  0  \ 

\pu  \  +  V  •  {pu)u  +  Vp  =  f  (1) 

\Ej^  \V-{Eu)J  \V-ipu)J 

where  we  have  split  the  flux  terms  into  an  advection  and  non-advection  part  and  lumped  viscous  terms  into 
f.  The  advection  part  (as  well  as  any  body  forces)  is  integrated  explicitly  to  give  intermediate  values  p*, 
{pHy  and  £’*.  Since  pressure  does  not  affect  the  continuity  equation,  =  p*.  The  momentum  update 

equation  can  be  divided  by  p^~^^  to  obtain 


(2) 

and  taking  its  divergence  gives 

(3) 

In  the  case  of  incompressible  flow,  we  would  set  V  •  =  0,  but  for  compressible  flow  we  instead  use  the 

pressure  evolution  equation  (see  e.g.  [9]), 


Pt  u  ’  Vp  =  —p(?V  •  u.  (4) 

If  we  fix  V  •  fx  to  be  at  time  through  the  time  step  (making  an  0{At)  error),  discretize  pt  +  u  •  Vp 

explicitly  using  a  forward  Euler  time  step  (i.e.  - — — h  •  Vp’^),  and  define  the  advected  pressure  as 

pa  —  pTi  _  obtain 


pU+l  =pa_  Atpc^^  •  U 

Substituting  this  in  Equation  (3)  and  rearranging  gives 


-n+l 


pn+l  _  -  /j”(c2)”AtV  •  V*, 


(5) 


(6) 


where  we  have  defined  pc^  at  time  and  the  pressure  p  at  time  Discretizing  the  gradient  and  divergence 
operators  yields 

1 


An+l 


G 


p 


n+l 


=  p^+p^(c")^AtG^r 


(7) 


where  G  is  our  discretized  gradient  operator,  —G^  is  our  discretized  divergence  operator,  and  p  and  u 
represent  variables  interpolated  to  cell  faces.  This  is  solved  to  obtain  p^+^  at  cell  centers.  The  time 
^n+i  pressures  are  then  applied  in  a  flux-based  manner  to  the  intermediate  momentum  and  energy  values  to 
obtain  time  quantities  in  a  discretely  conservative  manner  (thereby  giving  correct  shock  speeds).  This 

is  done  by  averaging  the  pressures  to  cell  faces  by  P^^i/2  ~  ^ ^ 


rewriting  Equation  (2)  using 

Pi  +Pi+i 
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face- averaged  quantities  Ui-^1/2  =  '^i+1/2  ~  At  —  (where  pi^ij2  =  {pi  +  Pi+i)/2),  and  updating  the 

values  using 


=  (pit)*  -  At 


^  n+l 
^i+1/2 


Pi-112 


Ax 


^nxi  -  At 


(p«)rfi/2-(p«)r-7/2 


Ax 


(8) 


3.  Solid  evolution 

We  give  a  brief  treatment  of  solid  evolution  with  sufficient  detail  to  properly  handle  the  fluid-structure 
interactions.  A  solid  state  is  completely  described  by  its  velocity  and  position.  We  update  the  position 
and  velocities  in  a  Newmark  scheme  in  which  velocity  at  time  is  used  to  update  the  position  to  time 

^n+i  ^  second  order  update.  Velocity  is  then  updated  from  time  t^  to  time  in  a  separate  step.  We 
describe  below  the  velocity  update  for  deformable  and  rigid  solids.  The  same  procedure  is  used  twice,  once 
with  a  time  step  of  At/2  to  obtain  for  position  update  and  then  with  a  time  step  of  At  for  the  final 

velocity  update. 

Deformable  body  formulation:  For  deformable  body  evolution  we  need  to  handle  both  elastic  and 
damping  forces.  Damping  forces  can  impose  strict  time  step  restrictions  and  are  thus  treated  implicitly.  We 
will  describe  a  method  which  treats  the  elastic  forces  explicitly  and  damping  forces  implicitly  although  one 
could  also  incorporate  implicit  elasticity.  The  deformable  body  at  a  given  time  t  can  be  described  by  a  vector 
of  positions  of  its  nodes  Xs{t)  and  a  vector  of  velocities  of  its  nodes  W(t).  The  evolution  of  velocities  can 
be  described  by  Newton’s  second  law  as 


Ms{Vs)t  =  F{Xs,Vs),  (9) 

where  Ms  is  the  mass  matrix  and  F  is  the  vector  of  all  forces  acting  on  the  solid  nodes.  Discretizing  and 
computing  the  elastic  terms  explicitly  and  damping  terms  explicit  in  position,  but  implicit  in  velocity,  i.e. 
F{X,,V,)  =  F{X^,Vr+^),  we  obtain 

+  AtF{X^,  (10) 

Using  a  Taylor  series  expansion  on  F  yields 

=  M,  W  +  At{F{X2,  V7)  +  -  y")).  (11) 

where  D  =  F{X^,Vp)  —  DV^  represents  the  elastic  only  (and,  if  present,  any  non-linear  damping 

terms  [19])  component  of  the  force  and  one  can  write 

=  M,V*  -f  (12) 

where  U/  denotes  the  velocity  vector  updated  explicitly  with  the  elastic  terms  only. 

Rigid  body  formulation:  For  a  rigid  body  we  define  the  generalized  velocity  vector  as  Vs  =  , 

where  Wm  is  the  velocity  of  its  center  of  mass  and  lo  is  its  angular  velocity.  The  velocity  evolution  can  then 
be  described  as 


where  is  a  3  x  3  diagonal  matrix  with  the  rigid  body  mass  in  the  diagonals,  Ir  is  the  inertia  tensor  and 
/,  r  are  the  net  force  and  torque  acting  on  it.  Writing  the  mass  matrix  as  Ms  and  combining  /,  r  into  F,  we 
get  a  form  similar  to  (9)  which  can  be  discretized  using  forward  Euler  to  obtain 

MsV^"^^  =  MsV^  +  AtF^  =  M,U/.  (14) 

Where  U/  denotes  the  velocity  vector  updated  with  the  explicit  forces.  Note  that  this  is  the  same  as 
Equation  (12)  except  without  any  damping  term.  We  will  therefore  use  Equation  (12)  as  our  general  solid 
update  equation  below,  as  it  covers  both  the  rigid  and  deformable  cases. 
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(a)  Eulerian  fluid  grid. 


(b)  Lagrangian  solid 
which  overlaps  the  fluid 
domain. 


(d)  Solid  nodes  which 
contribute  to  the  raster¬ 
ized  face. 


1_ 


(c)  Solid  voxelized  to 
fluid  faces. 


Figure  1:  A  common  challenge  with  FSI  problems  is  one  of  overlapping  grids.  We  resolve  this  issue  by  voxelizing  solid 
degrees  of  freedom  to  the  fluid  grid  using  an  interpolation  operator  denoted  by  the  matrix  W.  The  row  corresponding 
to  a  fluid  face  gets  contributions  from  nearby  solid  nodes. 


4.  Fluid-structure  interaction 

We  solve  for  the  fluid  on  an  Eulerian  grid,  and  the  solids  on  freely  deforming  Lagrangian  meshes.  The 
fluid  and  structure  interact  with  each  other  by  applying  equal  and  opposite  forces  at  the  interface,  satisfying 
physical  boundary  conditions  (we  use  no-slip,  no  penetration  boundary  conditions)  in  the  process.  Immersed 
boundary  methods  induce  extra  force  variables  at  the  interface  and  apply  a  regularization  operator  to  map 
these  forces  to  fluid  faces  (see  e.g.  [20]).  They  also  incorporate  an  interpolation  operator  to  map  fluid 
velocity  to  solid  nodes  for  applying  boundary  conditions.  We  eliminate  the  extra  interface  force  variables 
and  conservatively  map  the  fluid  pressures  directly  to  solid  nodes,  and  solid  velocities  to  fluid  faces  using  an 
interpolation  operator. 

Figure  1  illustrates  an  example  fluid  grid  which  is  coupled  to  a  Lagrangian  solid  which  occupies  the 
upper  right-hand  corner  of  the  grid.  In  our  model,  the  fluid  interacts  with  a  voxelized  version  of  the  solid 
and  the  solid  directly  sees  forces  acting  on  its  nodes.  We  define  an  interpolation  operator  W  which  maps 
solid  node  velocity  to  the  fluid  cell  faces,  where  the  rows  correspond  to  fluid  faces  and  the  columns  to  solid 
nodes.  W  can  be  constructed  in  a  row-by-row  fashion:  for  each  row,  we  identify  the  corresponding  fluid  face 
and  locate  the  nearby  solid  nodes.  The  entry  corresponding  to  each  solid  node  is  populated  by  a  weight 
proportional  to  its  contribution  to  the  fluid  face,  and  then  finally  the  row  is  normalized  to  ensure  that  each 
row  sums  to  one,  making  it  an  interpolation.  This  is  done  in  a  component-by-component  manner,  e.g.  the 
x-component  of  solid  velocity  is  voxelized  to  x-axis  fluid  faces  but  not  y-  or  2:-axis  fluid  faces,  and  so  the 
solid  velocity  at  fluid  face  i  +  1/2  is  (IEEs)^+i/2.  Since  pressure  is  defined  at  cell  centers,  we  also  introduce 
an  extrapolation  operator  B  which  maps  cell-centered  pressure  to  face  pressures,  as  illustrated  in  Figure  2. 
These  face  pressures  are  then  multiplied  by  the  surface  area  of  the  cell  face  to  get  a  force  and  distributed 
back  to  solid  nodes  using  .  That  is,  W  maps  from  solid  node  degrees  of  freedom  to  cell  faces,  and 
maps  back  in  the  opposite  direction.  Note  that  since  the  rows  of  W  sum  to  one,  the  columns  of  sum  to 
one  and  therefore  the  force  felt  due  to  the  pressure  on  the  face  is  fully  and  conservatively  distributed  to  the 
solid  node  degrees  of  freedom. 

4.1.  The  strongly  coupled  system 

The  fluid  acts  on  solid  degrees  of  freedom  via  pressure  along  the  interface.  The  pressure  exerts  a  force 
given  by  A f  Bp  on  the  solid  degrees  of  freedom,  where  Aj  is  a  diagonal  matrix  whose  entries  correspond 

to  the  areas  of  fluid- structure  faces.  We  can  incorporate  these  forces  into  the  implicit  solid  system  given  by 
Equation  (12): 

M,  =  MsV*  +  AtDV^+^  +  AtW'^A  fBp.  (15) 

The  fluid  sees  a  velocity  boundary  condition  at  the  fluid-structure  interface.  To  incorporate  this  into 
the  fluid  equations,  we  partition  the  discrete  divergence  operator  —G^  into  two  components.  Gj  operates 
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B  = 


0  0 
0  0 
0  1 


Figure  2:  Operator  B  maps  pressure  from  cell  centers  to  bordering  fluid-structure  faces.  In  this  example  there  are 
a:-direction  faces,  of  which  the  one  to  the  far  right  represents  a  rasterized  solid  face.  Therefore  B  has  three  rows 
(one  for  each  vertical  face,  with  the  top  and  the  bottom  rows  corresponding  to  the  far  left  and  far  right  vertical  faces 
respectively,  and  the  middle  row  corresponding  to  the  middle  vertical  face),  and  two  columns  (one  for  each  pressure 
at  each  cell  center).  Since  the  only  contribution  to  the  solid  is  from  the  second  pressure  to  the  third  face,  B  has  the 
form  shown  above  with  a  single  non-zero  element.  Note  that  {l/dx)B^  equals  —G^,  as  defined  in  Figure  3(b). 


over  fluid-fluid  faces,  while  Gj  is  the  component  of  the  divergence  operator  which  operates  on  rasterized 
fluid-structure  faces  (as  outlined  in  Figure  3),  and  =  Gj  +  Gj.  We  can  then  set  fluid-structure  faces  to 
have  implicit  Neumann  boundary  conditions;  that  is, 


u*  - 

P 


at  a  fluid- fluid  face;  and 
at  a  fluid-structure  face. 


Taking  the  divergence  of  the  velocity  field  yields 

^  qT ^yn+1 


(16) 


(17) 


Using  this  modified  definition  for  in  Equation  (5)  and  substituting  into  Equation  (3)  gives 


- — j  +  -JG 

Ai/9”(c")2  J 


pn+i_Gjwvr^  =  ^  +  G^u^ 


(18) 


If  we  define  V  =  AxAyAz  to  be  the  volume  of  the  fluid  cell,  then  VGj  =  AfB^.  Combining  equations 
(15)  and  (18),  using  scaled  pressure  p  =  Atp  and  scaled  advected  pressure  p^  =  Atp",  and  rescaling  the  fluid 
equations  by  cell  volume  gives  us  our  symmetric  system 


^I^VGJlGf  -AfB^W  \  ^  (^P^  +  VGju^\ 

-W^BAf  -Ms  +  AtDj\yj^^^J  \  -MsV^  )' 


(19) 


It  is  interesting  to  note  that  if  we  take  the  incompressibility  assumption  (i.e.  c  ^  oo)  then  this  system 
reduces  to  one  similar  to  [17]. 

The  system  in  Equation  (19)  is  symmetric  but  indefinite,  and  can  be  solved  using  efficient  solvers  such 
as  Conjugate  Residuals  [13]  to  obtain  the  final  time  solid  velocity  and  pressure.  The  solid  part  of  our 
update  is  now  complete,  but  we  still  need  to  use  the  pressure  to  update  the  fluid  momentum  and  energy 

(noting  that  is  already  done). 


4.2.  Updating  fluid  momentum  and  energy 

To  obtain  correct  shock  speeds  we  use  the  flux-based  method  discussed  above,  with  modifications  to 
account  for  fluid-structure  faces.  At  a  fluid-structure  face  i  +  1/2,  the  fluid  applied  a  force  of  {BAfp)i_^i/2  to 
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-1/dx  1/dx 

-1/dx 

(a)  Gj 


j_ 

dx 


-1 

0 


1  0 

-1  0 


(b) 


j_ 

dx 


0 

0 


0 
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Figure  3:  In  our  derivation,  the  divergence  operator  —G^  is  split  into  (which  operates  only  on  fluid- fluid  faces) 
and  G^  (which  operates  only  on  fluid-structure  faces).  We  show  this  splitting  for  a  simple  two  cell  example  where 
the  right-most  face  is  a  fluid-structure  interface.  The  rows  in  the  above  matrices  correspond  to  cells  and  columns 
to  faces.  The  left  most  face  corresponds  to  the  first  column  of  G^  and  only  has  one  non-zero  element  since  it  only 
borders  one  fluid  cell.  The  middle  face  (which  corresponds  to  the  second  column  of  Gj)  contributes  to  both  fluid 
cells  and  hence  has  two  non-zero  elements.  The  third  column  of  G^  is  zero,  as  the  third  face  is  a  fluid- structure  face 
and  instead  corresponds  to  GJ.  Figure  (b)  depicts  G'^ ,  which  is  dehned  as  —{l/dx)B^  in  Figure  2. 


the  solid.  To  conserve  momentum,  fluid  face  i  +  1/2  should  apply  an  equal  and  opposite  force  — (5^/p)i+i/2 
on  fluid  cell  i.  In  our  momentum  update  this  is  numerically  equivalent  to  setting  Pi^i/2  =  (^p)i+i/2 
fluid- structure  faces. 

Next,  we  need  to  consider  the  work  done  by  the  fluid  on  the  solid  at  a  fluid-structure  face.  We  are 
applying  an  impulse  on  the  solid,  which  is  equivalent  to  applying  a  constant  force  over  the 

interval  At.  In  order  to  compute  the  work  done  on  the  solid  system  by  a  single  force  /  in  the  presence  of 
other  forces,  we  lump  all  forces  acting  on  the  solid  into  a  vector  F  and  examine 


pAt 

Jo 


pAt 

Jo 


f  .  Vs{t)dt  =  /  •  (K"  +  M-^Ft)dt  =  Atf  . 


M7^fA 


=  Atf- 


V"  +  W+i 


(20) 


where  we  take  advantage  of  F  and  /  being  constant  over  the  interval.  We  are  interested  in  calculating  the 

work  done  by  a  single  fluid  face  on  the  solid,  so  if  we  take  kiT^i/2  column  vector  which  distributes 

the  pressure  from  cell  face  i  +  1/2  to  the  solid  node  degrees  of  freedom  then  /  =  W^^i2{BAfp)i_^ij2^  and 

the  work  done  on  the  solid  by  this  face  is  exactly 


At 


^2+1/2 


1  ^ 

- 

2 

—  At  [{BAfp)i_^i/2]  Wi_^i/2 


vr  + 


(21) 


This,  if  Pi+i/2  is  defined  to  be  (Bp) i_^i^2  suggested  above  in  the  momentum  update,  then  we  merely  need 
to  set  Ui^i/2  =  (l/2)(kk[l^/^  +  ill  order  to  obtain  a  flux  pu  which  exactly  conserves  the  kinetic 

energy  transferred. 


4.3.  Time  step  restrietion 

In  our  method  fluid-structure  interactions  are  handled  implicitly  and  thus  we  avoid  introducing  any  new 
time  step  restrictions.  The  time  step  is  therefore  determined  by  the  minimum  of  the  time  steps  imposed  by 
the  fluid  and  the  structure.  For  the  structure  update  the  time  step  restriction  is  determined  by  the  elastic 
part  only,  as  damping  terms  are  handled  implicitly,  while  our  semi-implicit  fluid  update  imposes  a  time  step 
restriction  dependent  only  on  its  bulk  velocity.  The  time  step  restriction  imposed  by  the  semi-implicit  flow 
formulation  in  two  spatial  dimensions  is 


l^lmaa;  \  ^  ^  |  ^  ^\Py\ 

Ay  J  pAx  pAy 


At  I  I'lf'lmaa:  f  \'^\max 

T  I  Ax  Ay  ^  V  V 
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and  we  refer  the  interested  reader  to  [11],  which  motivates  this  formulation. 

We  note  that  the  implicit  fluid- structure  coupling  gives  stable  results  even  for  very  high  density-to-mass 
ratios,  where  explicit  methods  struggle  even  when  the  CFL  restrictions  of  both  solid  and  fluid  systems  are 
obeyed.  We  explore  this  in  example  6.1.1. 


5.  Unified  time  integration 

We  employ  a  time  integration  scheme  which  incorporates  fluid  evolution  into  a  Newmark-style  solid 
evolution  scheme.  The  scheme  works  by  computing  an  intermediate  velocity  for  the  solid  and 

applying  this  in  a  second  order  update  to  get  solid  positions  at  time  .  Velocities  are  then  updated  from 
time  to  (discarding  intermediate  values),  and  so  two  linear  systems  are  solved. 

In  order  to  compute  the  intermediate  solid  velocity  we  begin  by  applying  all  explicit  solid  forces 

to  the  system,  which  gives  Explicit  body  forces  such  as  gravity  and  viscosity  are  also  applied  to 

the  fluid  system,  yielding  fluid  quantities.  The  coupled  system  (19)  is  solved  in  order  to  obtain 

j^n+i  _  ^  and  then  the  entire  fluid  state  and  all  solid  velocities  are  restored  to  their  time 

values. 

These  new  positions  are  then  used  to  compute  an  effective  velocity  for  the  solids,  i.e.  —  Xf)/ At. 

Using  the  effective  velocity  and  then  the  time  position  of  the  solid,  we  fill  ghost  cells.  These  ghost  cells  are 
used  directly  in  the  stencils  of  high-order  methods,  and  provide  a  valid  state  for  which  to  populate  uncovered 
cells.  In  order  to  compute  the  ghost  cell  data  at  location  Xg^  we  begin  by  identifying  the  closest  solid  interface 
point  x/,  and  reflecting  across  the  interface.  Density  and  pressure  are  interpolated  to  the  reflected  point 
2xi  —  Xg  from  neighboring  cells  and  then  copied  to  the  ghost  cell.  The  surface  normal  N  at  the  interface  is 
used  to  decompose  the  velocity  at  the  reflected  point  Vr  into  its  normal  component  VrN  =  Vr  •  N  and  its 
tangential  component  Wt  =  W  ~  WatZV.  In  order  to  remain  continuous  with  the  effective  velocity  of  the 
structure  at  the  interface  U/,  VrN  is  reflected  across  the  interface,  and  so  we  compute  VgN  =  2U/  '  N  —  VrN- 
Tangential  velocity  is  decoupled  from  the  interface  and  thus  we  can  use  it  directly,  giving  the  final  ghost  cell 
velocity  Vg  =  VgNX  +  VrT- 

Once  ghost  cells  are  filled,  explicit  body  forces  such  as  gravity  and  viscosity  are  integrated  into  the  system, 
and  the  advection  component  of  flux  from  Equation  (1)  is  applied  using  a  conservative  flux-based  method 
(see  [11].  Explicit  solid  forces  are  applied  in  order  to  compute  and  then  the  coupled  system  (19) 

is  solved  to  obtain  Vf^^  and  .  This  pressure  is  applied  as  per  Section  4.2  to  obtain  time  fluid 
quantities. 

We  also  fill  the  ghost  cells  inside  the  solid  using  time  data  from  the  fluid  and  solid  velocities,  as 
described  above.  Although  none  of  our  examples  use  these  ghost  values,  if  an  explicit  body  force  such  as 
viscosity  were  to  be  applied,  its  stencil  would  require  valid  ghost  cells  to  be  defined.  Note  that  these  are 
valid  as  instantaneous  ghost  cells,  whereas  the  ghost  cells  above  use  the  effective  solid  velocity,  which  is  the 
actual  motion  of  the  solid  through  the  mesh.  Practical  experience  shows  that  this  can  make  a  meaningful 
difference. 


6.  Examples  and  validation 

In  order  to  compare  our  results  with  previous  methods,  we  implement  an  explicit  coupling  scheme  which 
integrates  a  fully  explicit  compressible  flow  evolution  with  a  Newmark  time  integration  for  solids.  This 
explicit  method  proceeds  in  a  fashion  similar  to  Section  5,  except  that  instead  of  solving  the  system  (19)  we 
simply  fill  ghost  cells  inside  the  solid  once  and  explicitly  evolve  the  fluid  once,  while  time  pressures  along 
the  fluid- structure  interface  are  applied  to  the  solid  as  explicit  forces.  This  gives  us  an  explicitly  coupled 
time  evolution  scheme,  such  as  the  one  described  in  [7]. 

Although  one  might  assume  that  the  implicit  solve  would  cause  efficiency  bottlenecks,  we  observed  rela¬ 
tively  few  Conjugate  Residuals  iterations  per  time  step.  This  is  likely  due  to  the  strongly  diagonally  dominant 
nature  of  Equation  (19),  and  the  good  initial  guess  for  pressure  provided  by  the  equation  of  state  at  time 
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YoY  all  of  our  one  dimensional  examples  the  maximum  number  of  iterations  required  per  time  step  was 
3.  For  the  two  dimensional  examples,  the  rigid  body  coupling  example  required  a  maximum  of  4  iterations, 
while  the  deformable  coupling  example  required  a  up  to  24  iterations  per  time  step. 

In  all  of  the  examples  we  consider  the  fluid  is  simulated  using  an  ideal  gas  law,  with  7  =  1.4. 


6.1.  One- dimensional  validation 

We  examine  several  one  dimensional  fluid-structure  interactions  to  validate  our  method.  A  third  order 
ENO  scheme  [18]  is  used  along  with  an  advection-based  CFL  number  of  .6.  All  quantities  below  are  in  SI 
units,  with  density  as  kg/w? .,  pressure  in  Pa,  lengths  in  m,  spring  coefficients  in  Njm,  etc. 


6.1.1.  Sod  shook  eoupled  with  a  rigid  body 

Our  first  example  is  a  Sod  shock  interacting  with  a  rigid  body,  with  open  boundary  conditions.  The 
initial  condition  for  the  fluid  is 


(/9(a;,0),u(a;,0),p(a;,0)) 


(1,0,1)  if:c<.5, 

(.125,0,.!)  if:c>.5. 


A  rigid  body  of  mass  1  and  width  .2  starts  at  rest  with  its  center  of  mass  a  distance  of  .8  from  the  left  of  the 
domain.  The  domain  is  of  length  2.  The  rigid  body  remains  at  rest  until  the  shock  hits  it,  at  which  point  it 
accelerates  by  virtue  of  the  pressure  difference.  The  solid  body  continues  to  accelerate  until  it  converges  to 
a  velocity  of  .927453,  which  is  precisely  the  interfacial  velocity  of  the  Sod  Riemann  problem.  Figure  4  shows 
snapshots  of  the  pressure  profile  at  various  times  through  the  simulation.  For  comparison,  results  with  the 
explicit  method  are  shown  in  Figure  5.  We  also  do  a  convergence  analysis  of  our  method  in  Figure  6.  The 
error  in  the  position  of  the  rigid  body  is  computed  at  time  .9  from  the  highest  resolution  grid  simulated, 
which  is  6401  grid  cells.  The  convergence  order  of  the  error  is  estimated  as  1.6. 

It  is  interesting  to  consider  this  simple  problem  for  a  variety  of  density-to-mass  ratios.  Figure  7(a)  shows 
the  velocity  of  the  rigid  body  as  a  function  of  time  for  a  range  of  rigid  body  masses  in  the  semi-implicit 
case.  Figure  7(b)  shows  this  in  the  explicit  case.  We  note  that  the  explicit  simulation  struggles  with  high 
density- mass  ratios.  In  particular  it  appears  as  though  the  rigid  body  gains  too  much  momentum  in  a  single 
time  step,  causing  the  fluid  on  the  other  side  to  over- compress,  leading  to  a  very  stiff  oscillatory  system,  even 
though  the  time  step  obeyed  CFL  restrictions.  We  show  snapshots  of  the  pressure  profile  of  simulations  with 
a  light  solid  of  mass  .0001,  with  semi-implicit  and  explicit  schemes  in  Figure  8  and  Figure  9,  respectively. 


6.1.2.  Sod  shook  interaeting  with  a  fluid  piston 

We  consider  a  similar  problem,  this  time  with  solid  wall  boundary  conditions  and  a  larger  domain,  with 
the  initial  discontinuity  located  at  distance  1  from  the  left  of  the  domain.  The  rigid  body  has  a  mass  of  1, 
width  .2  and  starts  at  rest  with  its  center  of  mass  at  1.5  from  the  left  of  the  domain.  The  domain  is  of  length 
3.  The  shock  imparts  momentum  to  the  rigid  body  which  in  turn  compresses  the  fluid  on  its  right.  This 
compressed  fluid  creates  a  high  pressure  region  which  pushes  back  on  the  solid,  in  effect  creating  a  “fluid 
spring.”  This  causes  the  rigid  body  to  oscillate  as  shown  in  Figure  10,  which  plots  the  position  of  the  center 
of  mass  of  the  rigid  solid  as  a  function  of  time.  Figure  11  shows  snapshots  of  the  pressure  profile  at  various 
times  through  the  simulation.  For  comparison,  results  with  the  explicit  method  are  shown  in  Figure  12.  We 
also  do  a  convergence  analysis  of  our  method  in  Figure  13.  The  error  in  the  position  of  the  rigid  body  is 
computed  at  time  45  from  the  highest  resolution  grid  simulated,  which  is  6401  grid  cells.  The  convergence 
order  of  the  error  is  estimated  as  1.03. 


6.1.3.  Sod  shook  eoupled  with  a  mass-spring  system 

To  conclude  the  one-dimensional  examples,  we  consider  the  mass-spring  system  interacting  with  a  high 
pressure  gas  described  in  [1]  in  order  to  provide  validation  for  our  approach  against  an  analytic  solution. 
The  domain  is  of  length  20,  and  a  spring  is  fixed  to  the  right  side  of  the  domain  which  has  a  rest  length  of 
1,  a  stiffness  of  10^,  no  damping  and  a  mass  of  3.  The  fluid  is  given  by 

{p,p,u)  =  i4, 10®,  0) 
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Figure  4:  Semi-implicit  simulation  of  a  Sod  shock  hitting  a  rigid  body  of  mass  1.  Pressure  profile  of  the  fluid  is  shown  at  various 
times  through  the  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line  segment  at  the  bottom  of  the  plot,  with  pressure 
inside  the  solid  shown  as  a  linear  pressure  profile.  The  simulation  was  done  on  a  grid  of  resolution  1601. 
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(a)  t  =  0 


(b)  t  =  .25 


(c)  t  =  .b 


id)  t  =  i 


Figure  5:  Explicit  simulation  of  a  Sod  shock  hitting  a  rigid  body  of  mass  1.  Pressure  profile  of  the  fluid  is  shown  at  various 
times  through  the  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line  segment  at  the  bottom  of  the  plot,  with  pressure 
inside  the  solid  shown  as  a  linear  pressure  profile.  The  simulation  was  done  on  a  grid  of  resolution  1601. 
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Figure  6:  Position  error  of  the  center  of  mass  of  a  rigid  body  hit  by  a  Sod  shock,  as  compared  to  a  high-resolution  simulation,  at 
time  .9s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the  resolution  of  the  underlying  grid.  The  convergence 
rate  is  1.6. 


An  outflow  boundary  condition  is  used  for  the  left  side  of  the  domain.  The  spring  starts  at  rest  length  and 
is  compressed  by  the  gas.  Figure  14  shows  snapshots  of  the  pressure  profile  at  various  times  through  the 
simulation.  The  position  of  the  moving  end  of  the  spring  as  a  function  of  time  is  shown  in  Figure  15(a),  and 
a  convergence  analysis  in  Figure  15(b).  The  error  in  the  position  of  the  free  end  of  the  spring  is  computed 
at  time  .008,  and  is  compared  against  the  analytic  solution  provided  in  [1].  The  convergence  order  of  the 
error  is  estimated  as  1.16. 

6.2.  Two-dimensional  validation 

In  this  section  we  validate  our  method  for  the  multidimensional  case,  and  briefly  describe  a  symmetric 
positive-definite  reformulation  of  the  Equation  (19).  We  consider  interactions  with  both  rigid  and  deformable 
solids.  A  second  order  ENO  scheme  was  used  along  with  an  advection-based  CEL  number  of  .6. 

6.2.1.  Rigid  Cylinder  lift-off 

This  example,  which  is  suggested  by  [5,  10,  1],  examines  the  interaction  of  a  Mach  3  shock  with  a  rigid 
cylinder  initially  at  rest  on  the  floor  of  a  rectangular  channel.  The  cylinder  is  lifted  by  the  shock,  due  to  an 
asymmetric  reflection  of  the  incident  wave.  The  test  domain  is  1  x  .2,  with  the  initial  shock  front  positioned 
at  .08  from  the  left  boundary  and  the  remaining  domain  is  filled  with  the  gas  at  pressure  1  and  density  1.4. 
The  top  and  bottom  of  the  domains  are  rigid  walls,  the  left  boundary  is  fixed  to  be  the  post  shock  state  and 
an  outflow  boundary  condition  is  used  for  the  right  boundary.  The  cylinder  has  a  density  of  10.77,  a  radius 
of  .05  and  is  initially  located  at  (.15,  .05).  Eigure  16  shows  the  snapshot  of  the  simulation  for  a  selection  of 
times.  Our  results  compare  favorably  to  those  shown  in  [1],  and  converges  at  a  rate  of  .93. 

6.2.2.  Deforming  ey Under  lift-off 

This  example  is  similar  to  the  one  described  above  (in  Section  6.2.1),  except  that  the  rigid  cylinder 
is  replaced  by  a  deformable  mass-spring  system  with  222  triangles,  and  edge-  and  altitude-springs  with  a 
stiffness  of  .3.  The  density  of  the  sphere  is  10.77,  has  a  radius  of  .05  and  the  center  of  mass  is  initially 
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(a)  Semi-Implicit. 


(b)  Explicit. 


(c)  Semi-Implicit  symmetric  positive-definite  formulation. 


Figure  7:  Velocity  of  a  I-D  rigid  body  hit  by  a  Sod  shock,  as  a  function  of  time.  Simulations  were  done  on  a  grid  of  resolution 
1601.  All  simulations  were  run  with  a  CFL  number  of  .6,  where  the  explicit  simulation  CFL  is  based  on  |n|  ±  c  and  the 
semi-implicit  simulation  was  run  with  the  CFL  condition  specified  in  Equation  (22).  The  explicit  simulations  grow  increasingly 
unstable  as  mass  tends  to  zero,  giving  unusable  results  when  mass  reaches  .0001  (these  results  are  shown  in  Figure  9),  and 
crashes  for  lighter  masses.  As  mass  tends  to  zero,  the  momentum  absorbed  by  the  solid  tends  to  zero  and  the  shock  passes 
through  the  solid  relatively  unperturbed,  and  so  the  fiat  line  to  which  solid  velocities  appear  to  converge  is  in  fact  the  post-shock 
velocity  of  the  fiuid. 
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(a)  t  =  0 


(b)  t  =  .25 


(c)  t  =  .5 


id)  t  =  i 


Figure  8:  Semi-implicit  simulation  of  a  Sod  shock  hitting  a  light  solid  of  mass  .0001.  Pressure  profile  of  the  fluid  is  shown  at 
various  times  through  the  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line  segment  at  the  bottom  of  the  plot.  The 
simulation  was  done  on  a  grid  of  resolution  1601.  For  this  light  mass,  the  post-shock  state  remains  practically  undisturbed  as 
very  little  momentum  transfers  to  the  solid. 
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(a)  t  =  0 


(b)  t  =  .25 


(c)  t  =  .5 


(d)  t  =  1 


Figure  9:  Explicit  simulation  of  a  Sod  shock  hitting  a  light  solid  of  mass  .0001.  Pressure  profile  of  the  fluid  is  shown  at  various 
times  through  the  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line  segment  at  the  bottom  of  the  plot.  The  simulation 
was  done  on  a  grid  of  resolution  1601.  The  CFL  number  for  this  simulation  is  .6,  and  we  use  the  standard  compressible  flow 
CFL,  based  on  |n|  ±  c.  Despite  satisfying  a  reasonable  CFL  time  step  restriction,  a  fully  explicit  simulation  generates  unstable 
results,  and  even  goes  unstable  and  crashes  for  masses  lighter  than  .0001. 
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Figure  10:  The  position  of  the  piston  (Section  6.1.2)  is  plotted  as  a  function  of  time. 


located  at  (.15,  .05).  Figure  17  shows  snapshots  of  the  simulation  for  a  selection  of  times.  As  the  shock 
front  passes  through  the  deforming  body,  it  dissipates,  scatters  and  is  partially  absorbed  by  the  body.  The 
example  converges  at  a  rate  of  .99. 

6.2.3.  Heavy  deforming  ey Under  lift-off 

We  next  consider  a  heavy  deforming  cylinder,  in  the  same  setup  as  described  in  Section  6.2.1  and  Sec¬ 
tion  6.2.2  above.  In  this  case,  the  cylinder  matches  the  cylinder  from  Section  6.2.2,  except  the  density  is 
set  to  100.  As  the  body  absorbs  the  shock  wave,  it  compresses  and  delays  the  shock.  Some  of  the  shock  is 
reflected,  but  most  of  the  shock  passes  through  the  cylinder.  Figure  18  shows  snapshots  of  the  simulation 
at  a  selection  of  times.  The  example  converges  at  a  rate  of  1.01. 

6.2.4-  Shoek  traveling  down  a  deformable  tube 

This  example  is  similar  to  the  inflatable  bladder  examples  suggested  in  [1]  and  [7]  in  which  a  shock  wave 
travels  through  a  deformable  tube  causing  large  deformation  of  the  walls.  Our  results  are  shown  in  Figure  19. 
We  also  do  a  convergence  analysis  of  our  method  in  Figure  20.  The  error  in  the  position  of  a  particle  on 
the  deformable  tube  is  computed  at  time  .000495  (which  is  the  approximate  time  of  maximum  deformation 
of  that  particle  in  the  highest  resolution  simulation)  from  the  highest  resolution  grid  simulated,  which  is 
800  X  600  grid  cells.  The  convergence  order  of  the  error  is  estimated  as  1.18. 

6.2.5.  Symmetrie  positive- definite  reformulation 

Our  numerical  method  is  symmetric,  but  not  positive-definite.  Recent  developments  in  [16]  discuss  a 
modification  of  the  implicit  coupling  methodology  for  incompressible  flow  by  separating  out  the  coupling 
forces  as  implicit  variables  A  (similar  to  immersed  boundary  methods),  decomposing  the  symmetric  damping 
force  into  D  =  C^C  and  solving  for  W  =  The  symmetric  positive-definite  system  they  obtain  can 

be  modified  for  compressible  flow  in  a  manner  similar  to  Section  4.1  to  obtain 
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(a)  t  =  0 


(b)  t  =  1.5 


(c)  t  =  2.87  (d)  t  =  4 


Figure  11:  Semi-implicit  simulation  of  a  piston  hit  by  a  Sod  shock,  with  closed- wall  boundary  conditions  on  both  sides.  Pressure 
profile  of  the  fluid  is  shown  at  various  times  through  the  semi-implicit  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line 
segment  at  the  bottom  of  the  plot,  with  pressure  inside  the  solid  shown  as  a  linear  pressure  profile.  The  simulation  was  done 
on  a  grid  of  resolution  1601.  The  shock  on  the  left  pushes  the  rigid  body  and  compresses  the  fluid  on  the  right  into  a  small 
high  pressure  pocket  against  the  wall,  which  in  turn  pushes  the  rigid  body  back  to  the  left. 
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(a)  t  =  0  (b)  t  =  1.5 


(c)  t  =  2.87  (d)  t  =  4 

Figure  12:  Explicit  simulation  of  a  piston  hit  by  a  Sod  shock,  with  closed-wall  boundary  conditions  on  both  sides.  Pressure 
profile  of  the  fluid  is  shown  at  various  times  through  the  explicit  simulation.  The  1-D  rigid  body  is  drawn  as  a  blue  line  segment 
at  the  bottom  of  the  plot,  with  pressure  inside  the  solid  shown  as  a  linear  pressure  profile.  The  simulation  was  done  on  a  grid 
of  resolution  1601.  The  shock  on  the  left  pushes  the  rigid  body  and  compresses  the  fluid  on  the  right  to  a  very  high  pressure 
against  the  wall,  which  in  turn  pushes  the  rigid  body  back  to  the  left. 
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Figure  13:  Position  error  of  the  center  of  mass  of  the  piston  (Section  6.1.2),  as  compared  to  a  high-resolution  simulation,  at 
time  4s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the  resolution  of  the  underlying  grid.  The  convergence 
rate  is  1.03. 


/  +  0  \Ip\  (  +  \ 

-Kf3-^G  +  WM-^W^)K^  KWM-^G^  A  =  KWV^  -  Ku^  ,  (23) 

\  0  GM-^W^K^  I  +  GM-^G^  J  \  J  \  ) 

where  G  and  —G^  are  the  volume  weighted  gradient  and  divergence  operators  respectively,  l3  is  the  diagonal 
matrix  of  fluid  dual  cell  masses,  and  is  the  matrix  of  Is  and  Os  mapping  A  to  the  appropriate  fluid 
velocity  scalars  (see  [16]  for  more  details).  Note  that  in  order  to  avoid  confusion  in  notation  we  renamed  a 
few  operators.  In  particular  W  and  J  in  [16]  correspond  to  the  K  and  W  we  use  here,  respectively.  This 
system  is  both  symmetric  and  positive-definite.  We  demonstrate  the  viability  of  this  modified  method  in 
another  example,  where  we’ve  replaced  the  implicit  coupled  solve  with  Equation  (23).  Our  example  is  similar 
to  the  example  in  Section  6.2.1  except  that  the  sphere  is  replaced  with  a  diamond  whose  major  axis  is  of 
length  .1  and  minor  axis  is  of  length  .025.  The  diamond  begins  rotated  by  7r/4,  with  a  center  of  mass  at 
(.15,  .04).  Snapshots  of  the  resulting  simulation  are  shown  in  Figure  21.  The  convergence  analysis  for  this 
example  is  shown  in  Figure  22  which  estimates  the  convergence  order  of  the  error  as  .84. 

7.  Conclusions  and  future  work 

We  have  presented  a  first  order  method  which  implicitly  couples  compressible  flow  with  solid  bodies  with 
arbitrary  constitutive  models.  We  show  that  this  method  is  robust,  numerically  conservative,  and  avoids 
the  numerical  instabilities  which  comparable  explicit  methods  suffer  from  in  the  presence  of  high  density- 
to-mass  ratios.  The  same  methodology  can  be  applied  to  reformulate  our  implicit  system  into  a  symmetric 
positive- definite  system. 

There  are  several  interesting  avenues  of  future  work  which  we  wish  to  explore.  Given  the  promising  results 
which  arise  from  handling  fluid-structure  interactions  implicitly,  we  believe  that  an  alternative  approach 
would  split  the  fluid  flux  along  Riemann  invariants-rather  than  by  pressure-and  solve  for  the  Riemann 
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(a)  t  =  0 


(b)  t  =  .0015 


(c)  t  =  .003  (d)  t  =  .0045 


(e)  t  =  .0075  (f)  t  =  .01 

Figure  14:  Semi-implicit  simulation  of  a  1-D  mass-spring  system  hit  by  a  Sod  shock  wave.  Pressure  profile  of  the  fluid  is  shown 
at  various  times  through  the  semi-implicit  simulation.  The  mass-spring  system  is  drawn  as  a  blue  line  segment  at  the  bottom 
of  the  plot.  The  simulation  was  done  on  a  grid  of  resolution  1601.  Note  the  formation  of  a  spontaneous  shock  wave. 
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invariant  which  interacts  with  the  solid  implicitly.  Our  method  also  relies  on  the  assumption  that  the  solid 
has  some  thickness  where  ghost  cells  can  be  filled,  and  we  believe  that  the  method  can  be  made  to  work  for 
thin  shell  structures  (such  as  parachutes).  Given  the  utility  of  the  scheme  proposed  in  [11]  in  handling  fluid- 
structure  interactions,  it  becomes  imperative  to  address  the  issues  of  that  original  scheme.  In  particular,  the 
implicit  component  of  the  method  is  overly  centrally- differenced,  which  tends  to  introduce  Gibbs  phenomena 
at  shocks.  It  would  be  better  to  add  upwind  biasing,  although  it  is  unclear  how  to  do  so. 
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19.16 


(a)  Position  of  the  free  end  of  the  spring,  as  a  function  of  time. 


(b)  Position  error  for  the  left-most  side  of  the  mass-spring  system,  as  compared  to  the  analytic 
solution  provided  in  [1],  at  time  .008s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the 
log  of  the  resolution  of  the  underlying  grid.  The  convergence  rate  is  1.16. 


Figure  15:  1-D  mass-spring  system  hit  by  a  Sod  shock  wave. 
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(d)  Position  error  of  the  center  of  mass  of  the  cylinder  hit  by  a  planar  shock, 
as  compared  to  a  high- resolution  simulation,  at  time  t  =  .15s,  with  a  con¬ 
vergence  of .96. 


Figure  16:  Pressure  contours  for  semi-implicit  simulation  of  rigid  cylinder  lift  off  are  shown  at  t  =  0,  t  =  .164  and  t  =  .301. 
The  simulation  is  run  with  a  CFL  number  of  .6,  using  the  CFL  restriction  discussed  in  Equation  22. 
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(d)  Position  error  of  the  center  of  mass  of  the  deformable  cylinder  hit  by  a 
planar  shock,  as  compared  to  a  high-resolution  simulation,  at  time  t  =  .15s. 
We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the  resolution 
of  the  underlying  grid.  The  convergence  rate  is  .99. 


Figure  17:  Pressure  contours  for  semi-implicit  simulation  of  deformable  cylinder  lift  off  are  shown  ai  t  =  0,  t  =  .164  and 
t  =  .301.  The  simulation  is  run  with  a  CFL  number  of  .6,  using  the  CFL  restriction  discussed  in  Equation  22. 
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(d)  Position  error  of  the  center  of  mass  of  the  heavy  deformable  cylinder 
hit  by  a  planar  shock,  as  compared  to  a  high-resolution  simulation,  at  time 
t  =  .15s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the 
resolution  of  the  underlying  grid.  The  convergence  rate  is  1.01. 


Figure  18:  Pressure  contours  for  semi-implicit  simulation  of  deformable  cylinder  lift  off  are  shown  at  t  =  0,  t  =  .164  and 
t  =  .301.  The  simulation  is  run  with  a  CFL  number  of  .6,  using  the  CFL  restriction  discussed  in  Equation  22. 


26 


:.V\V: 

j'r  ^  'j : :::::::::::::  i: :::::::  :j  : 

i;  '5^ r r r r r ::::  ::  r r ;;;  r r r r . 
r  H  ^^3  J'?!  in  r  3  3 1 ;  3  j  j  r  r  3 !  1  r  3 !  1  r  •'  3 ;  r  j  3 ; ;  ? ;  1 ;  3 ; ;  r  ;  r  i ;  r  r  3 !  -  ' ' ' '  . . 

j  ^  w . 

if  ;  U  r  r  3 : :  1 : : :  r  r : ;  ^  r " ; : : : :  i :  r ; : :  r : :  r  r  n :  r : : :  r  3  3 :  r ;  i  ^  r  r : :  ^  r : : : : : : : : :  j : : 

f j  1  If : :  r: r; :::  r':: :  — 

f  1 1  1  i  5  s :  i :  L  ■  L  ■ : :  ^  ■  L  “  :  L : :  1  :  1 : 1 L  ■  :  L "  ■  “ : :  ■  “ ;  ■  “  :  1 1  ‘  J  M  L  ■: : :  i  : 

■  E J  :3iiJ : lll:: : l: :: F.: : L3 :: if: :  14.: :  1'^ ; 

1 f:::  r  I  ^ F. : ; 

=f",-  ■  -  ■ 

Figure  19:  A  planar  shock  travels  down  a  deformable  bladder.  Shown  are  the  velocity  field  of  the  fiuid  in  green  and  the  velocities 
of  the  deformable  nodes  in  red  at  times  t  =  .0001,  t  =  .0002,  t  =  .0003,  t  =  .0004,  t  =  .0005  and  t  =  .0006. 
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Figure  20:  Position  error  of  the  position  of  a  particle  on  the  deformable  tube  hit  by  a  planar  shock,  as  compared  to  a  high- 
resolution  simulation,  at  time  .00049s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the  resolution  of  the 
underlying  grid.  The  convergence  rate  is  1.18. 
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Figure  21:  A  diamond  is  hit  by  a  planar  shock,  and  then  collides  with  the  top  of  the  channel.  Shown  are  pressure  contours  at 
t  =  0^  t  =  .04,  t  =  .08,  t  =  .16  and  t  =  .2. 
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Figure  22:  Position  error  of  the  center  of  mass  of  the  diamond  hit  by  a  planar  shock,  as  compared  to  a  high-resolution 
simulation,  at  time  .15s.  We  plot  the  log  of  the  relative  error,  as  a  function  of  the  log  of  the  resolution  of  the  underlying  grid. 
The  convergence  rate  is  .84. 
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